准备数据

library(animalcules)
library(SummarizedExperiment)
library(MultiAssayExperiment)
#--得到phyloserq对象并提取必要数据信息
library(ggClusterNet)
library(phyloseq)
data(ps)

otu = as.data.frame(t(vegan_otu(ps)))
head(otu)
##           KO1  KO2 KO3  KO4  KO5  KO6  OE1  OE2  OE3 OE4  OE5  OE6  WT1  WT2
## ASV_1    1113 1968 816 1372 1062 1087 1270 1637 1368 962 1247 1017 2345 2538
## ASV_1591   12   10   1   14    3   12   15   16    4  16   31   18   17   58
## ASV_657     5    7   6    9    8    6   16   26    9   6    9   11    2    1
## ASV_28    253  303  54  439  132  182  154  263  153 257  242  178  391  450
## ASV_1717    3    6   0    1    0    4    2    6    2   0    1    2    4    0
## ASV_2407    1    2   7    0    1    7    2    2    1   2    3    8    4    1
##           WT3  WT4  WT5  WT6
## ASV_1    1722 2004 1439 1558
## ASV_1591   17   34   99   11
## ASV_657     3    8   18   14
## ASV_28    235  321  449  242
## ASV_1717    1    1    0    8
## ASV_2407    2    0    1    3
tax = as.data.frame((vegan_tax(ps)))
head(tax)
##           Kingdom         Phylum          Class           Order
## ASV_1    Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_1591 Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_657  Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_28   Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_1717 Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_2407 Bacteria Actinobacteria Actinobacteria Actinomycetales
##                       Family          Genus    Species
## ASV_1    Thermomonosporaceae     Unassigned Unassigned
## ASV_1591 Thermomonosporaceae     Unassigned Unassigned
## ASV_657  Thermomonosporaceae     Unassigned Unassigned
## ASV_28   Thermomonosporaceae Actinocorallia Unassigned
## ASV_1717 Thermomonosporaceae   Actinomadura Unassigned
## ASV_2407 Thermomonosporaceae   Actinomadura Unassigned
map = sample_data(ps)
head(map)
##     Group      Date    Site Sequencing  Platform     Species Batch
## KO1    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO2    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO3    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO4    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO5    KO  2017/7/4  Harbin        BGI HiSeq2500 Arabidopsis     1
## KO6    KO  2017/7/5  Harbin        BGI HiSeq2500 Arabidopsis     1
##     BarcodeSequence LinkerPrimerSequence      ReversePrimer
## KO1      ACGCTCGACA  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO2      ATCAGACACG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO3      ATATCGCGAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO4      CACGAGACAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO5      CTCGCGTGTC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO6      TAGTATCAGC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
#--首先构造SummarizedExperiment对象,比较简单,类似phyloseq对象
micro <- SummarizedExperiment(assays=list(counts=as.matrix(otu)),
                              colData=map,
                              rowData=tax)
# 将SummarizedExperiment对象封装成为ExperimentList

mlist <- ExperimentList()
mlist[[1]] = micro
names(mlist) = "MicrobeGenetics"# 注意必须命名,否则无法区分每个部分数据组
# 构造不同数据组之间的记录文件
gistmap <- data.frame(
  primary = row.names(map),
  colname = row.names(map),
  stringsAsFactors = FALSE)
maplistowe <- list(MicrobeGenetics = gistmap)
sampMapowe <- listToMap(maplistowe)

# colData文件为分组文件,数据框即可,本案例只有一个微生物组数据,所以直接用map文件就可以了。
#-下面就直接构建了MultiAssayExperiment文件
mae <- MultiAssayExperiment(experiments = mlist, colData = map,
                            sampleMap = sampMapowe)

运行Shiny app-导入数据-网页分析

微生物组数据基本统计和描述

这部分用于统计每个样本中OTU的数量,并做两种方式可视化:频率曲线,箱线 + 散点图;如果使用shiny程序的话,直接可以展示表格。 此外,可以按照微生物分类水平合并OTU数据:

p <- filter_summary_pie_box(mae,
                            samples_discard = c("subject_2", "subject_4"),
                            filter_type = "By Metadata",
                            sample_condition = "Group")
p

更换分组,重新统计。

p <- filter_summary_bar_density(mae,
                                samples_discard = c("subject_2", "subject_4"),
                                filter_type = "By Metadata",
                                sample_condition = "Group")
p

微生物组数据分析-微生物丰度展示-堆叠柱状图

通过tax_level选择某个分类等级,通过 sample_conditions 选择需要添加的分组标签。值得注意的是这里可以对堆叠柱状图排序的,通过order_organisms来指定,默认丰度从高到低。这里从源代码来看就是通过改变factor来实现的,所以图例第一个也是排序的这个微生物。

p <- relabu_barplot(mae,
                    tax_level="Family",
                    # order_organisms=c('Retroviridae'),
                    sort_by="nosort",
                    sample_conditions= "Group",
                    show_legend=TRUE)
p

微生物组数据分析-微生物丰度展示-热图

p <- relabu_heatmap(mae,
                   tax_level="Genus",
                   # sort_by="conditions",
                   sample_conditions=c("Group"))
p

微生物丰度展示-箱线图

library(tidyverse)
tax = vegan_tax(ps) %>% as.data.frame()
Gname <- table(tax$Genus) %>%names()
p <- relabu_boxplot(mae,
                    tax_level="Genus",
                    organisms=Gname[1:3],
                    condition="Group",
                    datatype="logcpm")
p

微生物组数据分析-多样性-alpha多样性

这里有四个多样性指标。然后通过箱线+散点图展示。

alpha_div_boxplot(MAE = mae,
                  tax_level = "Genus",
                  condition = "Group",
                  alpha_metric = "shannon") 

对多样性进行统计检验。这里可选的是“Wilcoxon rank sum test”, “T-test”, “Kruskal-Wallis”这三种方法。这里作者检测仅仅支持分组文件是数值的(因为作者使用了cor.test函数),不够人性化。

otu = as.data.frame(t(vegan_otu(ps)))
head(otu)
##           KO1  KO2 KO3  KO4  KO5  KO6  OE1  OE2  OE3 OE4  OE5  OE6  WT1  WT2
## ASV_1    1113 1968 816 1372 1062 1087 1270 1637 1368 962 1247 1017 2345 2538
## ASV_1591   12   10   1   14    3   12   15   16    4  16   31   18   17   58
## ASV_657     5    7   6    9    8    6   16   26    9   6    9   11    2    1
## ASV_28    253  303  54  439  132  182  154  263  153 257  242  178  391  450
## ASV_1717    3    6   0    1    0    4    2    6    2   0    1    2    4    0
## ASV_2407    1    2   7    0    1    7    2    2    1   2    3    8    4    1
##           WT3  WT4  WT5  WT6
## ASV_1    1722 2004 1439 1558
## ASV_1591   17   34   99   11
## ASV_657     3    8   18   14
## ASV_28    235  321  449  242
## ASV_1717    1    1    0    8
## ASV_2407    2    0    1    3
tax = as.data.frame((vegan_tax(ps)))
head(tax)
##           Kingdom         Phylum          Class           Order
## ASV_1    Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_1591 Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_657  Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_28   Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_1717 Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_2407 Bacteria Actinobacteria Actinobacteria Actinomycetales
##                       Family          Genus    Species
## ASV_1    Thermomonosporaceae     Unassigned Unassigned
## ASV_1591 Thermomonosporaceae     Unassigned Unassigned
## ASV_657  Thermomonosporaceae     Unassigned Unassigned
## ASV_28   Thermomonosporaceae Actinocorallia Unassigned
## ASV_1717 Thermomonosporaceae   Actinomadura Unassigned
## ASV_2407 Thermomonosporaceae   Actinomadura Unassigned
map = sample_data(ps)
head(map)
##     Group      Date    Site Sequencing  Platform     Species Batch
## KO1    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO2    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO3    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO4    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO5    KO  2017/7/4  Harbin        BGI HiSeq2500 Arabidopsis     1
## KO6    KO  2017/7/5  Harbin        BGI HiSeq2500 Arabidopsis     1
##     BarcodeSequence LinkerPrimerSequence      ReversePrimer
## KO1      ACGCTCGACA  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO2      ATCAGACACG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO3      ATATCGCGAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO4      CACGAGACAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO5      CTCGCGTGTC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO6      TAGTATCAGC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
map$Group = map$Group %>% as.numeric()

#--首先构造SummarizedExperiment对象,比较简单,类似phyloseq对象
micro <- SummarizedExperiment(assays=list(counts=as.matrix(otu)),
                              colData=map,
                              rowData=tax)
# 将SummarizedExperiment对象封装成为ExperimentList

mlist <- ExperimentList()
mlist[[1]] = micro
names(mlist) = "MicrobeGenetics"# 注意必须命名,否则无法区分每个部分数据组
# 构造不同数据组之间的记录文件
gistmap <- data.frame(
  primary = row.names(map),
  colname = row.names(map),
  stringsAsFactors = FALSE)
maplistowe <- list(MicrobeGenetics = gistmap)
sampMapowe <- listToMap(maplistowe)

# colData文件为分组文件,数据框即可,本案例只有一个微生物组数据,所以直接用map文件就可以了。
#-下面就直接构建了MultiAssayExperiment文件
mae <- MultiAssayExperiment(experiments = mlist, colData = map,
                            sampleMap = sampMapowe)



do_alpha_div_test(MAE = mae,
                  tax_level = "Genus",
                  condition = "Group",
                  alpha_metric = "inverse_simpson",
                  alpha_stat = "T-test")
##                                       output
## Method  Pearson's product-moment correlation
## P-value                     0.71148699581376

微生物组数据分析-多样性-beta多样性-聚类距离热图

diversity_beta_heatmap(MAE = mae, 
                       tax_level = 'Genus', 
                       input_beta_method = "bray",
                       input_bdhm_select_conditions = 'Group',
                       input_bdhm_sort_by = 'condition')

其次通过组内距离和组件距离的箱线图展示

diversity_beta_boxplot(MAE = mae, 
                       tax_level = 'Genus', 
                       input_beta_method = "bray",
                       input_select_beta_condition = 'Group')

再有就是统计检验,共有三种方法可以选择:PERMANOVA,Kruskal-Wallis,Wilcoxon test。 但是只有两种距离可供选择,其次就是两两比较不能实现。

diversity_beta_test(MAE = mae, 
                    tax_level = 'Genus',
                    input_beta_method = "bray",
                    input_select_beta_condition =  'Group',
                    input_select_beta_stat_method = 'PERMANOVA',
                    input_num_permutation_permanova = 999)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = dist.mat ~ condition, data = sam_table, permutations = input_num_permutation_permanova)
##           Df SumOfSqs      R2     F Pr(>F)  
## condition  1 0.029817 0.12357 2.256  0.019 *
## Residual  16 0.211468 0.87643               
## Total     17 0.241285 1.00000               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

微生物组数据分析-排序 四种方法

PCA

result <- dimred_pca(mae,
                     tax_level="Genus",
                     color="Group",
                     shape="Group",
                     pcx=1,
                     pcy=2,
                     datatype="logcpm")
result$plot

PCoA

result <- dimred_pcoa(mae,
                      tax_level="Genus",
                      color="Group",
                     shape="Group",
                      axx=1,
                      axy=2,
                      method="bray")
result$plot

UMAP

esult <- dimred_umap(mae,
                      tax_level="Genus",
                      color="Group",
                     shape="Group",
                      cx=1,
                      cy=2,
                      n_neighbors=15,
                      metric="euclidean",
                      datatype="logcpm")
result$plot

t-SNE,除了二维图形展示还可以进行三维图形的展示。

result <- dimred_tsne(mae,
                      tax_level="Genus",
                      color="Group",
                     shape="Group",
                      k="3D",
                      initial_dims=30,
                      perplexity=10,
                      datatype="logcpm")
result$plot

微生物组数据分析-差异分析

p <- differential_abundance(mae,
                            tax_level="Phylum",
                            input_da_condition=c("Group"),
                            min_num_filter = 2,
                            input_da_padj_cutoff = 0.5)
p
##                        microbe     padj   pValue log2FoldChange experiment
## 1               Proteobacteria 5.66e-05 3.54e-06          1.030     1: 6/6
## 2               Planctomycetes 2.08e-04 2.60e-05         -1.350     1: 6/6
## 3               Actinobacteria 4.05e-03 9.18e-04          0.591     1: 6/6
## 4              Verrucomicrobia 4.05e-03 1.01e-03         -0.852     1: 6/6
## 5                Acidobacteria 1.43e-01 5.35e-02         -0.552     1: 6/6
## 6                   Unassigned 1.43e-01 5.03e-02          0.480     1: 6/6
## 7              Armatimonadetes 1.49e-01 6.50e-02         -0.843     1: 6/6
## 8  Candidatus_Saccharibacteria 1.88e-01 9.38e-02          0.946     1: 6/6
## 9                  Chloroflexi 2.89e-01 1.63e-01          0.459     1: 6/6
## 10              Thaumarchaeota 4.65e-01 2.91e-01         -1.300     1: 2/6
## 11                  Firmicutes 3.40e-06 2.13e-07          1.570     1: 6/6
## 12              Proteobacteria 1.61e-02 2.01e-03          0.689     1: 6/6
## 13 Candidatus_Saccharibacteria 6.15e-02 1.48e-02          1.400     1: 6/6
## 14             Verrucomicrobia 6.15e-02 1.54e-02         -0.631     1: 6/6
## 15             Armatimonadetes 7.15e-02 2.23e-02         -1.040     1: 6/6
## 16                  Unassigned 1.08e-01 4.07e-02          0.502     1: 6/6
## 17              Planctomycetes 3.17e-01 1.39e-01         -0.482     1: 6/6
## 18                 Nitrospirae 4.26e-01 2.13e-01         -0.472     1: 6/6
## 19             Ignavibacteriae 4.42e-01 2.49e-01          0.988     1: 5/6
## 20              Actinobacteria 4.44e-01 2.80e-01          0.193     1: 6/6
## 21              Thaumarchaeota 4.44e-01 3.05e-01         -1.260     1: 2/6
## 22                  Chlamydiae 4.52e-01 3.39e-01         -0.453     1: 6/6
## 23               Acidobacteria 4.94e-01 4.32e-01         -0.225     1: 6/6
## 24                Spirochaetes 4.94e-01 4.12e-01         -0.287     1: 6/6
## 25                  Firmicutes 3.75e-07 2.34e-08          1.690    2: 12/6
## 26              Planctomycetes 4.53e-02 5.66e-03          0.872    2: 12/6
## 27              Actinobacteria 1.35e-01 2.54e-02         -0.399    2: 12/6
## 28                 Chloroflexi 3.89e-01 1.03e-01         -0.536    2: 12/6
## 29              Proteobacteria 3.89e-01 1.22e-01         -0.345    2: 12/6
## 30               Acidobacteria 4.43e-01 2.49e-01          0.327    2: 12/6
## 31             Ignavibacteriae 4.43e-01 2.18e-01          1.040    2: 11/6
## 32                 Nitrospirae 4.43e-01 1.93e-01         -0.483    2: 12/6
## 33                Spirochaetes 4.43e-01 2.46e-01         -0.401    2: 12/6
##    control prevalence Group Size adjusted fold change Contrast
## 1   2: 6/6     66.67%                            1.00  1 vs. 2
## 2   2: 6/6     66.67%                            1.00  1 vs. 2
## 3   2: 6/6     66.67%                            1.00  1 vs. 2
## 4   2: 6/6     66.67%                            1.00  1 vs. 2
## 5   2: 6/6     66.67%                            1.00  1 vs. 2
## 6   2: 6/6     66.67%                            1.00  1 vs. 2
## 7   2: 6/6     66.67%                            1.00  1 vs. 2
## 8   2: 6/6     66.67%                            1.00  1 vs. 2
## 9   2: 6/6     66.67%                            1.00  1 vs. 2
## 10  2: 5/6     38.89%                            2.50  1 vs. 2
## 11  3: 6/6     66.67%                            1.00  1 vs. 3
## 12  3: 6/6     66.67%                            1.00  1 vs. 3
## 13  3: 6/6     66.67%                            1.00  1 vs. 3
## 14  3: 6/6     66.67%                            1.00  1 vs. 3
## 15  3: 6/6     66.67%                            1.00  1 vs. 3
## 16  3: 6/6     66.67%                            1.00  1 vs. 3
## 17  3: 6/6     66.67%                            1.00  1 vs. 3
## 18  3: 6/6     66.67%                            1.00  1 vs. 3
## 19  3: 6/6     61.11%                            1.20  1 vs. 3
## 20  3: 6/6     66.67%                            1.00  1 vs. 3
## 21  3: 5/6     38.89%                            2.50  1 vs. 3
## 22  3: 6/6     66.67%                            1.00  1 vs. 3
## 23  3: 6/6     66.67%                            1.00  1 vs. 3
## 24  3: 6/6     66.67%                            1.00  1 vs. 3
## 25  3: 6/6    100.00%                            2.00  2 vs. 3
## 26  3: 6/6    100.00%                            2.00  2 vs. 3
## 27  3: 6/6    100.00%                            2.00  2 vs. 3
## 28  3: 6/6    100.00%                            2.00  2 vs. 3
## 29  3: 6/6    100.00%                            2.00  2 vs. 3
## 30  3: 6/6    100.00%                            2.00  2 vs. 3
## 31  3: 6/6     94.44%                            1.83  2 vs. 3
## 32  3: 6/6    100.00%                            2.00  2 vs. 3
## 33  3: 6/6    100.00%                            2.00  2 vs. 3

微生物组数据分析-生物标记物分析

这里可选的方法有两个:“logistic regression”, “random forest”。这里去除一个分组。因为ROC只能支持两个分组。

ps <- subset_samples(ps,!Group %in% c("OE"));ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2861 taxa and 12 samples ]
## sample_data() Sample Data:       [ 12 samples by 10 sample variables ]
## tax_table()   Taxonomy Table:    [ 2861 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2861 tips and 2860 internal nodes ]
otu = as.data.frame(t(vegan_otu(ps)))
head(otu)
##           KO1  KO2 KO3  KO4  KO5  KO6  WT1  WT2  WT3  WT4  WT5  WT6
## ASV_1    1113 1968 816 1372 1062 1087 2345 2538 1722 2004 1439 1558
## ASV_1591   12   10   1   14    3   12   17   58   17   34   99   11
## ASV_657     5    7   6    9    8    6    2    1    3    8   18   14
## ASV_28    253  303  54  439  132  182  391  450  235  321  449  242
## ASV_1717    3    6   0    1    0    4    4    0    1    1    0    8
## ASV_2407    1    2   7    0    1    7    4    1    2    0    1    3
tax = as.data.frame((vegan_tax(ps)))
head(tax)
##           Kingdom         Phylum          Class           Order
## ASV_1    Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_1591 Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_657  Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_28   Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_1717 Bacteria Actinobacteria Actinobacteria Actinomycetales
## ASV_2407 Bacteria Actinobacteria Actinobacteria Actinomycetales
##                       Family          Genus    Species
## ASV_1    Thermomonosporaceae     Unassigned Unassigned
## ASV_1591 Thermomonosporaceae     Unassigned Unassigned
## ASV_657  Thermomonosporaceae     Unassigned Unassigned
## ASV_28   Thermomonosporaceae Actinocorallia Unassigned
## ASV_1717 Thermomonosporaceae   Actinomadura Unassigned
## ASV_2407 Thermomonosporaceae   Actinomadura Unassigned
map = sample_data(ps)
head(map)
##     Group      Date    Site Sequencing  Platform     Species Batch
## KO1    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO2    KO 2017/6/30 Beijing        BGI HiSeq2500 Arabidopsis     1
## KO3    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO4    KO  2017/7/2   Sanya        BGI HiSeq2500 Arabidopsis     1
## KO5    KO  2017/7/4  Harbin        BGI HiSeq2500 Arabidopsis     1
## KO6    KO  2017/7/5  Harbin        BGI HiSeq2500 Arabidopsis     1
##     BarcodeSequence LinkerPrimerSequence      ReversePrimer
## KO1      ACGCTCGACA  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO2      ATCAGACACG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO3      ATATCGCGAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO4      CACGAGACAG  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO5      CTCGCGTGTC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
## KO6      TAGTATCAGC  AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC
#--首先构造SummarizedExperiment对象,比较简单,类似phyloseq对象
micro <- SummarizedExperiment(assays=list(counts=as.matrix(otu)),
                              colData=map,
                              rowData=tax)
# 将SummarizedExperiment对象封装成为ExperimentList

mlist <- ExperimentList()
mlist[[1]] = micro
names(mlist) = "MicrobeGenetics"# 注意必须命名,否则无法区分每个部分数据组
# 构造不同数据组之间的记录文件
gistmap <- data.frame(
  primary = row.names(map),
  colname = row.names(map),
  stringsAsFactors = FALSE)
maplistowe <- list(MicrobeGenetics = gistmap)
sampMapowe <- listToMap(maplistowe)

# colData文件为分组文件,数据框即可,本案例只有一个微生物组数据,所以直接用map文件就可以了。
#-下面就直接构建了MultiAssayExperiment文件
mae <- MultiAssayExperiment(experiments = mlist, colData = map,
                            sampleMap = sampMapowe)


p <- find_biomarker(mae,
                    tax_level = "Genus",
                    input_select_target_biomarker = c("Group"),
                    nfolds = 3,
                    nrepeats = 6,
                    seed = 99,
                    percent_top_biomarker = 0.2,
                    model_name = "logistic regression")
#-提取生物标记物
p$biomarker
##         biomarker_list
## 1          Thermogutta
## 2        Bacteriovorax
## 3         Ohtaekwangia
## 4    Methyloversatilis
## 5        Paenibacillus
## 6          Kineosporia
## 7         Chitinophaga
## 8           Lysobacter
## 9      Virgisporangium
## 10     Ammoniibacillus
## 11         Rhizobacter
## 12        Methylovorus
## 13                Gp11
## 14                 Gp6
## 15           Kribbella
## 16 Armatimonadetes_gp5
## 17      Permianibacter
## 18        Chryseolinea
## 19                Gp17
## 20        Aquihabitans
## 21           Rhizobium
## 22         Nannocystis
## 23       Brevundimonas
## 24          Haliangium
## 25         Turneriella
## 26           Niastella
## 27        Streptomyces
## 28  Methyloceanibacter
## 29            Opitutus
## 30            Bacillus
## 31      Bradyrhizobium
## 32           Pelomonas
## 33             Devosia
## 34        Sphingomonas
## 35            Massilia
## 36        Nocardioides
## 37       Piscinibacter
## 38      Steroidobacter
## 39      Hyphomicrobium
## 40         Polaromonas
## 41       Aeromicrobium
## 42         Ramlibacter
## 43    Phenylobacterium
#对重要变量可视化。
# importance plot
p$importance_plot

#ROC曲线准确度评估。注意ROC曲线只能对二分便量进行操作。
# ROC plot
p$roc_plot

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.936 
## [2] LC_CTYPE=Chinese (Simplified)_China.936   
## [3] LC_MONETARY=Chinese (Simplified)_China.936
## [4] LC_NUMERIC=C                              
## [5] LC_TIME=Chinese (Simplified)_China.936    
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] caret_6.0-88                lattice_0.20-44            
##  [3] forcats_0.5.1               stringr_1.4.0              
##  [5] dplyr_1.0.6                 purrr_0.3.4                
##  [7] readr_1.4.0                 tidyr_1.1.3                
##  [9] tibble_3.1.2                ggplot2_3.3.3              
## [11] tidyverse_1.3.1             phyloseq_1.36.0            
## [13] ggClusterNet_0.1.0          MultiAssayExperiment_1.18.0
## [15] SummarizedExperiment_1.22.0 Biobase_2.52.0             
## [17] GenomicRanges_1.44.0        GenomeInfoDb_1.28.0        
## [19] IRanges_2.26.0              S4Vectors_0.30.0           
## [21] BiocGenerics_0.38.0         MatrixGenerics_1.4.0       
## [23] matrixStats_0.58.0          animalcules_1.9.1          
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1           backports_1.2.1        igraph_1.2.6          
##   [4] plyr_1.8.6             lazyeval_0.2.2         splines_4.1.0         
##   [7] crosstalk_1.1.1        BiocParallel_1.26.0    digest_0.6.27         
##  [10] foreach_1.5.1          htmltools_0.5.1.1      fansi_0.5.0           
##  [13] magrittr_2.0.1         memoise_2.0.0          covr_3.5.1            
##  [16] cluster_2.1.2          limma_3.48.0           recipes_0.1.16        
##  [19] Biostrings_2.60.0      annotate_1.70.0        modelr_0.1.8          
##  [22] gower_0.2.2            askpass_1.1            prettyunits_1.1.1     
##  [25] colorspace_2.0-1       rvest_1.0.0            blob_1.2.1            
##  [28] haven_2.4.1            xfun_0.23              crayon_1.4.1          
##  [31] RCurl_1.98-1.3         jsonlite_1.7.2         genefilter_1.74.0     
##  [34] survival_3.2-11        iterators_1.0.13       ape_5.5               
##  [37] glue_1.4.2             gtable_0.3.0           ipred_0.9-11          
##  [40] zlibbioc_1.38.0        XVector_0.32.0         DelayedArray_0.18.0   
##  [43] Rhdf5lib_1.14.0        shape_1.4.6            GUniFrac_1.2          
##  [46] rentrez_1.2.3          scales_1.1.1           DBI_1.1.1             
##  [49] Rcpp_1.0.6             plotROC_2.2.1          progress_1.2.2        
##  [52] viridisLite_0.4.0      xtable_1.8-4           reticulate_1.20       
##  [55] bit_4.0.4              tsne_0.1-3             lava_1.6.9            
##  [58] umap_0.2.7.0           prodlim_2019.11.13     DT_0.18               
##  [61] glmnet_4.1-1           htmlwidgets_1.5.3      rex_1.2.0             
##  [64] httr_1.4.2             RColorBrewer_1.1-2     ellipsis_0.3.2        
##  [67] farver_2.1.0           pkgconfig_2.0.3        XML_3.99-0.6          
##  [70] dbplyr_2.1.1           nnet_7.3-16            sass_0.4.0            
##  [73] locfit_1.5-9.4         utf8_1.2.1             labeling_0.4.2        
##  [76] tidyselect_1.1.1       rlang_0.4.11           reshape2_1.4.4        
##  [79] AnnotationDbi_1.54.0   cellranger_1.1.0       munsell_0.5.0         
##  [82] tools_4.1.0            cachem_1.0.5           cli_2.5.0             
##  [85] generics_0.1.0         RSQLite_2.2.7          ade4_1.7-16           
##  [88] broom_0.7.6            evaluate_0.14          biomformat_1.20.0     
##  [91] fastmap_1.1.0          yaml_2.2.1             fs_1.5.0              
##  [94] ModelMetrics_1.2.2.2   knitr_1.33             bit64_4.0.5           
##  [97] KEGGREST_1.32.0        nlme_3.1-152           xml2_1.3.2            
## [100] rstudioapi_0.13        compiler_4.1.0         plotly_4.9.3          
## [103] png_0.1-7              reprex_2.0.0           geneplotter_1.70.0    
## [106] bslib_0.2.5.1          stringi_1.6.2          highr_0.9             
## [109] RSpectra_0.16-0        Matrix_1.3-3           multtest_2.48.0       
## [112] shinyjs_2.0.0          vegan_2.5-7            permute_0.9-5         
## [115] vctrs_0.3.8            pillar_1.6.1           lifecycle_1.0.0       
## [118] rhdf5filters_1.4.0     jquerylib_0.1.4        data.table_1.14.0     
## [121] bitops_1.0-7           R6_2.5.0               codetools_0.2-18      
## [124] MASS_7.3-54            assertthat_0.2.1       rhdf5_2.36.0          
## [127] openssl_1.4.4          DESeq2_1.32.0          withr_2.4.2           
## [130] GenomeInfoDbData_1.2.6 hms_1.1.0              mgcv_1.8-35           
## [133] grid_4.1.0             rpart_4.1-15           timeDate_3043.102     
## [136] class_7.3-19           rmarkdown_2.8          pROC_1.17.0.1         
## [139] lubridate_1.7.10

reference